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Abstract 

In this paper we investigate the acceleration of the expansion of pre- 
mixed spherical flames and evolution of the cellular patterns on their sur- 
faces. An asymptotic model is used for the simulations and a spectral numer- 
ical algorithm is employed to study flames over large time intervals. Numer- 
ous numerical experiments indicate that for large enough time intervals the 
acceleration of two-dimensional expanding flame ceases and the expansion 
rate stabilizes to a value significantly exceeding the burning rate. The impor- 
tance of the effect of forcing was also confirmed and the validity of sectorial 
simulations of closed flames was studied in order to justify prospective use 
of the Fourier spectral model for the three-dimensional spherical flames. 

1 Introduction 

In [1] extensive experimental data on premixed expanding spherical flames have 
been analysed and it was concluded that starting from a certain moment the av- 
eraged flame radius < r > grows as < r > oc t 3 ^ 2 rather than < r > oc t. The 
phenomenon was linked to the cellularization of flames, which was well known 
from experiments too. Indeed, the appearance of cellular patterns increases the 
flame surface area, hence the fuel consumption, and hence the averaged flame 
expansion rate. 

The cellularization of flame fronts was, in its turn, associated with intrinsic 
combustion instabilities. The effect of the hydrodynamic combustion instability 
on expanding spherical flames was studied in [2] using the linear perturbation 
theory combined with phenomenological assumptions. Later, the approach was 
further improved and freed from the phenomenological assumptions [3]. These 
linearized solutions confirmed the onset of the instability of the flame front but 
could not quantify its cellularization and acceleration because the latter phenom- 
ena are essentially nonlinear. 
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A simple, yet physically reasonable, nonlinear model of hydrodynamically 
unstable planar flames was suggested in [4]. In [5] it was extended to expanding 
spherical flames and studied numerically confirming that there is a time instance 
£*, such that the flame expansion rate behaves like <r> oc t 3 / 2 for t > £*. The 
same model was also suggested in [6], where similar to the planar flames [7], 
analytical pole solutions were obtained and studied as well. 

Because of their physical origins simple nonlinear models of expanding flames 
[5, 6] are expected to be valid only locally. Apparently, the results obtained when 
applying them to the whole flame are instructive indeed, but still inconclusive and 
cannot be accepted as the adequate theoretical model of cellularization and accel- 
eration. A physically consistent global model of flames of any smooth enough 
geometry was developed in [8]. Mathematically, the approach projects the gov- 
erning equations to the flame surface reducing the mathematical dimension of the 
problem by one. However, the resulting equation is still extremely costly from 
the computational point of view and only two-dimensional simulations have been 
carried out so far. 

A compromise between universality and computability was suggested in [9], 
where consideration was limited to a narrow but still very practical case of flames 
which do not deviate from the spherical ones significantly. On the technical side 
the model combines the operator of the linearized problem obtained in [2] for 
the expanding spherical flame in terms of spherical harmonics expansions and a 
Huygens type nonlinearity specific to the local nonlinear model [5, 6]. Physically, 
model [9] is consistent with [8] and is robust and plausible enough to simulate 
the cellularization of expanding spherical flames in three spatial dimensions. At 
the time of writing of this paper, the flame sizes we were able to reach in our 
computations do not significantly exceed those reported in [9] and are not large 
enough to match our two-dimensional calculations. However, our investigations 
show that numerical studies of the expanding three-dimensional flames on the 
time scales required for comparison with the two-dimensional calculations are 
possible. 

In the following sections we specify the mathematical models and numerical 
algorithms to solve them. In Section |4] we report our results on flame front be- 
haviour on long time intervals and on the effect of external forcing. Also, we 
present our attempts to simulate the three-dimensional flames and assess the pos- 
sibility of simulation of closed flames via their finite segments. 

2 Mathematical Models 

Let us consider an expanding flame front and assume that its surface is close 
enough to a sphere and that every point on the flame surface is uniquely defined 
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by its distance r = r(0, <p, t) from the origin for < < 7r, < (f) < 2n, and 
t > 0. It is convenient to represent such a flame as a perturbation $(0, 0, t) of a 
spherical surface of a reference radius r (t), i.e. r(9,(f),t) = r (t) + $(0,0, £). 
Then, the Fourier image of the governing equation of the flame front evolution in 
the nondimensional notations suggested in [5, 6] can be written as 
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Here |&| < oo, t > 0, are the Fourier components of the properly scaled 
upstream perturbations of the unburnt gas velocity field /(0, t), and initial values 
of $fc(0) = are given. By construction, equation ([TJ holds in the sector 
< 4> < 2ir/9 w with a large enough integer 6 n . 

Models of [5] and [6] differ by an additive term proportional to 5o,fc/r (t). 
This term adds just 0(m t) to for k = only. This is not essential and the term 
is not included in our conception. 

Equation dU) was obtained as a local model of a curved expanding flame. How- 
ever, we will use it globally, on the whole flame surface, with 9 W = 1. In order 
to justify such an action, let us first note that if transformed back to the physical 
space, the equation takes the form 
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which is rotation-invariant for /(</>, £) = and 0^ = 1. 

On the other hand, we may rewrite the geometrically invariant equation ob 
tained in [8] in the coordinate form as follows: 
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( [r(0)-r(0)]-n(0) 



n(0) 



\r((j)) -r(6)\ 2 

[r(0)] 2 + 2(r ) 2 -r 

2{[r(0)] 2 + (r > ) 2 } ]U 

[— (r^ sin + r cos 0) , [r^ cos — r sin . 
\A" 2 + i. r <t>) 2 



and 



>0, 



(4) 



(5) 
(6) 



r(0) = [r(0) cos0, r(0) sin0]. 

Assuming again that r(0, t) = r (t) + $(0, t) and linearizing all but the Huygens 
terms, one arrives to the equation 
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which differs from ©, for 9^ = 1, by a few non-essential terms only. 

Comparison of © and © reveals the detailed relationship between the sim- 
plified model [5,6] and the comprehensive one [8]. In particular, it validates use 
of as a global model of the whole spherical flame as long as the perturbations 
remain small enough. 

Equations of the three-dimensional model [9] can be written in terms of the 
spherical harmonics expansion coefficients 
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Here \n\, \m\ < oo, t > 0, are the spherical harmonics coefficients of 

the properly scaled upstream perturbations of the unburnt gas velocity field, and 
initial values of $ n ,m(0) = are given. The expression for the linear response 

2n(n-l) (rc-l)(n + 2) 
" (n ' t)= (2 n+ l)r„ W ' <9> 

emerges from the analysis of [2] and N n ,m{t) are the spherical harmonics coeffi- 
cients of the nonlinear Huygens term 

Choice of r (t) in the models introduced in this Section may be based on a 
variety of principles. In equation Q, which governs small perturbations of ©- 
©, choice of ro(i) is more or less arbitrary. However 
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ro(t) = — / r(0,t)# (11) 
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is probably the most appropriate one because it minimizes the perturbations r{t) — 
ro(t) = $(0, t). On the other hand, equation ©, or 0, was obtained as that one, 
which governs perturbations of a uniform and steadily expanding spherical flame 
r(t) = t. Therefore, choice of 

ro(t) = t (12) 

is more reasonable in this case, as this is the solution around which the lineariza- 
tion is performed. Equation (QJ was already considered with (fTTTl in [5]. Here we 
complement the analysis of © by considering it with (fT2b . 



3 Computational Algorithms 

System dU is solved numerically by neglecting the harmonics of orders higher 
than a finite integer number K > 0. Then, the nonlinearity can be represented 
as a circular convolution and evaluated effectively with the FFT Also, we found 
that the stability of the numerical integration scheme suggested in [5] can be im- 
proved significantly by calculating the contribution from the linear terms in ([TJ 
analytically. Thus, the linear terms, i.e. the source of physical instability, are 
tackled exactly and only the nonlinear ones, with the dumping effect, are approx- 
imated. This improvement allowed us to continue the calculations for up to ten 
times further in time than in [5]. 
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Using the notation 

2 K 

9M = ~fJ2 l ( k - l)*&k-i - V5 k , + f k (t), (13) 

l=-K 

equation CQ) can be written as 
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and we will search for its solutions in the form 
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where Y(t) is the solution of the uniform equation. 
Straightforward evaluations yield 
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where t is an arbitrary real of which the final result (1131) does not depend. 

In order to integrate the ODE's (fT71) over the interval [t n ,t n + At] or [t n — 
At, t n + At], we interpolate values of g k (t) = <7fc[$] polynomially: 
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For example, the first order extrapolation gives 

g k (t) = G ktnfi = g k (t n ), (19) 
and the second one results in 
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In what follows we provide formulas for the first order extrapolation only. For- 
mulas of higher orders of accuracy are a bit bulky, but can be obtained straightfor- 
wardly. 

The interpolation results, in general, in 
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In order to move any further, we will now assume that 

r (t) oc t a . 

Then (l23l can be written as Au k (t) = uj k (t) — cu k (t n ), where 
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Then, the integration formula 422b is transformed into 
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with Wfc(t) given by and (ETl 

A computational algorithm based on numerical spherical harmonics transfor- 
mation (SHT) [10], similar to [9], was used in this work in order to solve (EJ). In 
addition, the stability of the numerical integration scheme was improved by eval- 
uating the contribution from the linear terms analytically, as shown above, and the 
code was parallelized in order to speed up the computations and to use larger data 
sets. 
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4 Computational results 



4.1 Saturation on the long time intervals 

Typical shapes of the flame fronts governed by ([TJ over large time intervals are 
illustrated in Fig. [lj The graph of [r((j>, t)— < r >]/< r > for t = 7.65 x 10 4 
shows that the wrinkle amplitudes are up to 10% of the averaged flame radius. 
The explicit forcing was not applied in this example. 




1.0 0.5 0.5 1.0 ' 0.1 0.05 0.05 0.1 



r($,t) xlO 5 (r-<r>)/<r> 



Figure 1: Evolution of a spherical flame governed by (HJ). Here values of r(0, t) 
are on the left and [r(<f>, t)— <r>]/ <r> for t = 7.65 x 10 4 are on the right. 
Positive values of the latter are in blue and negative ones are in red; 7 = 0.8, 

f(M = o. 

Permanent growth of the size of a spherical flame as it expands prompts studies 
of the effect of the size of a planar flame on its propagation speed as the first 
step towards the understanding of the acceleration mechanism of the expanding 
flames. The investigation of the dynamics of planar flames shows that the spatially 
averaged flame speed 

L 

1 f d$ f 

<$ < >= L J ~df dX (30) 



of a flame size L stops growing and begins to oscillate irregularly around its time 
average for large enough t. A definite correlation between the size of the 

flame and its stabilized spatially averaged propagation speed was estab- 

lished, see e.g. [11]. The effect was explained by proving the high sensitivity of 
planar cellular flames to particular types of linear perturbations, see [11-14]. By 
continuing calculations reported in [11] for even larger planar flames, we estab- 
lished that their propagation speed no longer grows after a certain critical flame 
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size is reached, see graph on the right of Fig. |2j In this paper we are interested in 
extending these findings for planar propagating flames to the expanding ones. In 
particular, we are studying the possibility of a stabilization of the expansion rate 
for large enough time intervals, when the flame size grows sufficiently large. 

Stabilized velocities <5C $ 4 y> of the planar flames and averaged velocities 
<$ t >= 7^ Jq W ^§t-d(j) of the spherical flames are compared in Fig. |2| Power law 
approximations (t — for the expansion rate of the spherical flame are also 
depicted there. Sudden increase of <$f > begins from t c « 2.2 x 10 3 . For the 
whole considered time interval [2.2 x 10 3 , 7.65 x 10 4 ] the optimal a « 0.34. For 
earlier times t E [2.2x 10 3 , 2. Ox 10 4 ], the best approximation is with a ~ 0.47, i.e. 
almost 1/2 as obtained in experiments. However, as time goes by, the expansion 
rate slows down and for t G [3.0 x 10 4 , 7.65 x 10 4 ] we got a « 0.23. 




Figure 2: Averaged flame front velocity for spherical (left) and planar (right) 
flames versus time t and flame size L respectively. Here 7 = 0.8 and /(</>, t) = 0. 
Markers on the right graph show the calculated cases. 

A tendency towards stabilization of <$ t > to a constant for t — > 00 is evident 
in the graph too. A change in morphology of the flame front at t s ~ 2.5 x 
10 4 is even more obvious in the graph of the variation of the perturbation of the 
averaged circular flame relatively to the radius of this averaged flame, i.e. of 
max [r(0, £)— <r>]/<r>, which is shown on the left of Fig. |3] 

The number of cells per unit length of the averaged flame contour is depicted 
on the right of Fig. |3] It stabilizes to a cell size of about 200, which is exactly 
the same as for large enough planar flames. Eventually, Fig. |4] illustrates the 
spectral distribution of energy of the perturbation $(0, t), which stabilizes by t s ~ 
2.5 x 10 4 too. 

Stabilization of the averaged perturbation gradient V$ = <9$/<90 can also be 
seen on the left of Fig. Besides the saturation of the gradient the graph illus- 
trates the importance of the number of Fourier modes K involved in the numerical 
simulation. The jumps in the graph correspond to the instances when we doubled 
K in order to match the continuously increasing size of the flame. The graph is 
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Figure 3: Ratio of the maximal amplitude of perturbation and of the averaged 
flame radius (left) and number of cells per unit length of the averaged flame con- 
tour (right). Here 7 = 0.8 and /(</>, t) = 0. 




Figure 4: Spectral distributions of energy of the perturbation $(0, t) for t = 
2.4928 x 10 4 (left) and t = 7.65 x 10 4 (right). 

getting smoother if K is adjusted more gradually and is kept much larger than the 
ratio of the critical wavelength A cr = 4 71/7 to the circumferential length of the 
flame. The graph on the right of Fig. |5] depicts the energy of the highest order 
Fourier harmonic involved in the simulation. 




xlO~ 



xlO 




Figure 5: Averaged gradient of the perturbation (left) and energy of the highest 
order harmonic involved in the simulation (right). 
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According to the data obtained in numerical simulations we may summarize 
the mechanism of flame front expansion in the Fourier space as follows. The 
energy of relatively long wave harmonics A > A cr permanently grows according to 
the Darrieus-Landau instability. This gained energy is transferred towards shorter 
wavelengths via the nonlinear effects and dissipates through the modes of short 
enough wavelengths A < A cr . As flame expands the number of short enough 
angular modes reduces and the overall dissipation rate may become insufficient 
to counterbalance the generation of energy due to the Darrieus-Landau instability. 
This results in an accumulation of energy in the short wavelength spectrum, see 
Fig. 13 (right), and leads to the blow-up of the numerical solution if the number K 
of used harmonics is not increased in time. This effect explains a slight elevation 
of the short wave tail of the spectral energy distribution on the right of Fig. |4] 
in comparison to the graph on the left. The latter one corresponds to the time 
moment soon after K was doubled, though the former one is just approaching the 
moment when K needs to be increased. 

In general, the process of stabilization of the expanding spherical front to a 
saturated state is very similar to the planar flame. The only distinctive difference 
is that the transitional period for the expanding flame is much longer. Using the 
realistic set of dimensional parameters from [5], e.g. planar flame speed relative to 
the burnt gases Ub = 0.5 m/ sec and thermal diffusivity D t h = 2.5 x 10~ 5 m 2 / 'sec, 
one may interpret our findings in dimensional terms as follows. First cusps begin 
to appear on the flame surface for <r> « 0.1 m; flame acceleration with the rate 
<r> oc i 3 / 2 starts for <r> « 0.4 m; acceleration rate begins to slow down for 
<r> « 5 m, and the acceleration ceases for <r> 20 m. Most of experiments 
summarized in [1] and reported in more recent works [15, 16] were carried out 
in enclosures and were affected by essential pressure rise and acoustics. Thus, 
there might be a considerable discrepancy between the characteristic flame radii 
just reported and those measured in realistic combustion experiments. Eventually, 
of course there is a chance that the stabilization of the expansion rate is just a 
consequence of the perturbative nature of the governing model used in this work. 

4.2 Effect of forcing 

A random point-wise set of perturbations uniformly distributed in time and in the 
Fourier space is a suitable model for both the computational round-off errors and 
a variety of perturbations of physical origins. In general, such a model would look 
like 

M(t) 

f(x,t) = ^ a m cos(£ w + <p m )6(t - t m ), (31) 

m=l 
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are non-correlated random sequences. It is assumed 

< t M (t) < t, < (p m < 2tt, and £ TO > 0, 



where a m , t m , £ m , and (p Tl 
that h < t 2 < ■ ■ ■ < t m < 
m = 1,2, ... , M(t). However, in practice we use only two harmonics with £ = 
= floor [7 <r>/2] and £ = + 1 weighted according to their closeness to the 
critical wavenumber £ cr = 7<r>/2: 



f(x,t) « /o[(l-£cr- + C*)cOs£*0+(£ cr -£*)cOs(£* + l] 
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5(t-t m ). (32) 



These two harmonics approximate the critical planar flame harmonic of the wave- 
length A cr = 47r/7, which was shown to contribute the most in the fastest growing 
perturbations in the Sivashinsky-type models of flame dynamics, see [17]. 

The sign of the noise amplitude f in (1321 was either plus or minus for every 
m with the probability 1/2. The delta function 5(t — t m ) was approximated by 
(7rr) _1 / 2 e~^~' m ^ l T with a small enough value of r. The impulse-like noise (l32l) 
is used here for the sake of simplicity. Some arguments towards its validity were 
suggested in [18]. Validation of models like (I3TT) . d32l and of the present numer- 
ical approach was carried out in [17]. More sophisticated and physically realistic 
models of temporal noise characteristics can be used with (HJ) as well. 

The importance of the effect of the amplitude of noise and frequency of ap- 
pearance of the impulse-like perturbations on the flame speed is illustrated in Fig. 
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Figure 6: Effect of forcing amplitude and frequency on the averaged flame expan- 
sion speed for <t m+1 — t m >= 50 (left) and <t m+1 — t m >= 10 (right). 



More details of the effect of noise are presented in Fig. [7] On the left we 
plotted graphs of the time instance at which the flame begins to accelerate ver- 
sus frequency of appearance of the impulse-like perturbations for a variety of the 
perturbation amplitude fo- Similar graphs of the nearly saturated flame expansion 
speed are given on the right. 
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Figure 7: Effect of forcing amplitude and frequency on the beginning of acceler- 
ation (left) and stabilized flame expansion speed (right). 



4.3 Three-dimensional numerical experiments 

In order to integrate (l8t- ([T0l) we used an algorithm which is very similar to that 
one developed in [9]. In addition, we parallelized the algorithm and used a few 
techniques to improve its stability and accuracy, see [19]. 

The basis functions in the spherical harmonics transformation (SHT) are the 
orthonormal eigenf unctions of the Laplace operator in spherical coordinates. How- 
ever, unlike the Discrete Fourier Transformation (DFT), the discrete SHT of a data 
set is its approximation rather than an equivalent representation. Namely, it links 
a full matrix with iV 2 elements representing a function on a regular spherical mesh 
of N 2 nodes in physical space and a triangular matrix of N 2 /2 spherical harmon- 
ics coefficients. Thus, every back/forward cycle results in loss of information. 
This loss is similar to the effect of a short wave filter. In our approach it is associ- 
ated with the nonlinear term only. 

Coordinate singularities at the sphere poles result in accumulation of approx- 
imation errors near the poles and weaken numerical stability. In order to prevent 
accumulation of approximation errors near the poles, we rotate the coordinate sys- 
tem around an axis in the equatorial plane by an angle uj from time to time. The 
coordinate transformation formulas for such an axis passing through the points 
4> = ±7r/2 are 

cos 9' = cos 9 cos uj + sin 9 sin uj cos <p 

cot 4>' = cot 4> cos uj — esc <p sin uj cot 9 

Their structure precludes use of the addition theorem for spherical polynomials to 
transform the coefficients <&„ )iri to the new coordinate system without the global 
SHT back to the physical space. On the other hand, the back/forward SHT entails 
the application of a short wave filter to the whole solution $ n>TO , rather than just to 
the nonlinear term N n m . Hence, these rotations should not be done too frequently. 
Also, it is useful to combine them with rotations around the axis passing through 
the poles. Implementation of the latter ones is trivial. 
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Message Passing Interface (MPI) parallelization paradigm was implemented 
to allow the computational work to be distributed to a number of processors. How- 
ever, the distributed data needs to be exchanged between these processors, which 
creates the biggest problem in modelling dynamics of large radius flames. Un- 
like the multidimensional DFT, discrete SHT does not possess a structure of a 
tensor product of one-dimensional transformations. Therefore, the data transmis- 
sion required by the truly distributed discrete SHT is much more sophisticated and 
intensive than just back/forward row/column transposition of the global solution 
array required by the DFT, see e.g. [20]. This results in a very fast growth of the 
communication overheads when the number of processors increases. 

Spherical computational algorithms based on parallel Legendre-Fourier trans- 
formations are much less efficient than the standard planar Fourier methods. How- 
ever, expanding spherical flames can be successfully simulated by the parallel SFT 
method. An example of the evolution of a random three-dimensional perturbation 
of a spherical flame is illustrated in Fig. [8] 




Figure 8: Evolution of a spherical flame; t ~ 40 (left) and t ~ 72 (right). Here 
7 = 0.764. 



4.4 Dynamics of curved flame segments 

Our numerical experiments showed that long time simulations of the three-dimensional 
expanding flames using discrete SHT might be possible on a parallel computer 
with large enough physically shared memory. However, there is an alternative 
approach based on simulations of the three-dimensional flame segments and fol- 
lowing extension to the whole surface by periodicity. In order to validate this idea 
we applied it to the two-dimensional flames first. 
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Results of two-dimensional numerical simulations of the dynamics of sectors 
9 £ [0, 27r/6' 7r ] of a spherical flame are illustrated in Fig. |9l One may see that 
the narrowing of the sector does not affect the flame expansion rate in absence of 
explicit forcing and that a correlation between this rate and the size of the segment 
becomes apparent as forcing strengthens. Similar observations were obtained for 
other parameters discussed earlier in this Section. 




Figure 9: Averaged flame segment expansion speed for / = (left) and f = 0.1, 
<t m+ i - t m >= 10 (right). 

Our calculations also show that there is a critical value of 9 W above which sec- 
torial simulations no longer represent the whole spherical flame. For example, we 
would not present the results of our simulations with 9 n = 16 as an approximation 
for the whole spherical flame. However, based on our two-dimensional results, it 
looks like simulations of the three-dimensional spherical flames using the Fourier, 
rather than the Fourier-Legendre, spectral model of the Sivashinsky type in the 
sector < 9, <fi < n, or even in < 9, <$ < w/2, are safe. 



5 Conclusions 

Long time interval simulations of a simplified model of the expanding spherical 
flames indicated that their expansion rate slows down as the flame size grows. The 
saturation of the planar flame propagation speed as their size grows was estab- 
lished too. Hence, a hypothesis of stabilization of the spherical flame expansion 
rate over a finite time interval is proposed. 

Further similarities with the propagating planar flames achieved by the ex- 
panding spherical flames on large time intervals were established in the studies 
of the effect of forcing, revealing a clear correlation between the strength of the 
forcing and the flame expansion rate. This supports the idea that the acceleration 
of both planar and expanding flames results from explicit and/or implicit forcing, 
which is always present both in experiments and calculations at least as noise of 
various physical origins. 
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In spite of many benefits of the simplified flame dynamics models, they have 
not been constructed to cope with significant perturbations of spherical flames. 
This casts a reasonable doubt in the possibility of extending the effects observed 
for the simplified model to realistic flames. In order to verify the hypothesis, 
numerical simulations of a more sophisticated model are required. In particular, 
model [8] is valid for flames of any geometry if thermal gas expansion due to 
combustion does not generate significant vorticity, i.e. for 7 -C 1. However, even 
in the coordinate form Q-© the governing equation of model [8] is extremely 
difficult to solve numerically, because, in contrast to ©, its nonlocal term has no 
convolution structure. 

Using parallelized spherical harmonics transformation, the evolution of a three- 
dimensional expanding spherical flame has been successfully simulated to a stage 
when wrinkles appear and form a well developed cellular structure. However, 
computational problems associated with the spherical harmonics transformation 
make it difficult to extend these calculations on time intervals which would match 
those attained in our two-dimensional simulations. On the other hand, it was also 
noticed that the simulations of the three-dimensional spherical flames using the 
Fourier, rather than the Fourier-Legendre, spectral model of the Sivashinsky type 
in the sector < 9, (ft < tt, or even in < 9, (ft < tt/2, are reasonable. 
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